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Flagella beating in complex fluids are significantly influenced by viscoelastic stresses. Relevant 
examples include the ciliary transport of respiratory airway mucus and the motion of spermatozoa 
in the mucus-filled female reproductive tract. We consider the simplest model of such propulsion 
and transport in a complex fluid, a waving sheet of small amplitude free to move in a polymeric 
fluid with a single relaxation time. We show that, compared to self-propulsion in a Newtonian fluid 
occurring at a velocity Un, the sheet swims (or transports fluid) with velocity 

U_ = 1 + DeV/r? 
U N 1 + De 2 ' 

where n 3 is the viscosity of the Newtonian solvent, r\ is the zero-shear-rate viscosity of the polymeric 
fluid, and De is the Deborah number for the wave motion, product of the wave frequency by the 
fluid relaxation time. Similar expressions are derived for the rate of work of the sheet and the 
mechanical efficiency of the motion. These results are shown to be independent of the particular 
nonlinear constitutive equations chosen for the fluid, and are valid for both waves of tangential and 
normal motion. The generalization to more than one relaxation time is also provided. In stark 
contrast with the Newtonian case, these calculations suggest that transport and locomotion in a 
non-Newtonian fluid can be conveniently tuned without having to modify the waving gait of the 
sheet but instead by passively modulating the material properties of the liquid. 

I. INTRODUCTION 

G.I. Taylor all but started a new field half a century ago when he showed that a two-dimensional waving sheet 
could self-propel in a viscous fluid in the absence of inertia Much progress toward understanding small-scale 
biological locomotion has been made since Taylor's contribution, and a number of classical reviews address the 
fundamental fluid mechanical issues associated with swimming microorganisms 0,0,0,0]. In this paper, motivated 
by the observation that a variety of internal biological fluids have complex polymeric microstructures and display 
non-Newtonian behavior, we revisit Taylor's original calculation in the case where the fluid is non-Newtonian. 

We start by reviewing two examples where complex fluids are involved in biological locomotion and transport, 
namely ciliary transport of respiratory mucus and sperm motion in the mucus-filled female reproductive tract. After 
summarizing the modeling approaches used in the past to quantify the influence of viscoelasticity on these examples, 
we motivate and introduce the model used in this paper. 

A. Biological background 

1 . Ciliary transport of respiratory mucus 

The first example of a viscoelastic biological fluid is the mucus present in our upper respiratory tract. It is a 
glycoprotein-based gel 0, 0] which is transported along the epithelium by large arrays of beating cilia, allowing 
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removal of foreign substances and a healthy respiratory system [8[ . Cilia are short flagella which produce fluid motion 
by means of a collaborative beating process (metachronal waves) and have long been recognized as crucial in many 
biological transport processes, in particular for respiratory and reproductive functions, as well as sensing mechanisms 
and propulsion 

In the upper respiratory tract, cilia are located in a thin layer of low- viscosity Newtonian fluid (serous layer) and 
penetrate partially in the overlying high- viscosity non-Newtonian mucous blanket. The rheology of respiratory airway 
mucus is critical to ensure normal transport conditions and many respiratory diseases such as asthma or cystic fibrosis 
lead to rheological changes 10( . Early rheological measurements of respiratory mucus reported a large shear viscosity, 
on the order of 1000 Pas, and a relaxation time scale A w 100 — 1000 s lCj, lll|. Later measurements showed a larger 
range of viscosities, between 10 and 1000 Pas 0, with significantly shorter relaxation times, A ~ 30 s [l3 |. 
strong shear-thining behavior [3| and elastic modulus on the order of 1 — 10 Pa [lj[ [ljj [if]]. Given that common 
cilia oscillate at frequencies of about / = uj/2-k 5-50 Hz [3j, the fluid mechanics of ciliated transport in respiratory 
mucus occurs therefore in the limit where the Deborah number, Dc = Acj, is much larger than one, and elastic effects 
are therefore predominant. 



2. Motion of spermatozoa in the female reproductive tract 



The second example of non-Newtonian biological fluid arises when we consider the locomotion of spermatozoa in the 
female reproductive tract. During their journey to the ovum, sperm cells are both actively swimming and passively 
transported by peristaltic-like flows (inside the uterus). Along the way, they encounter a variety of viscoelastic fluids, 
through the cervix (glycoprotein-based cervical mucus), inside the uterine cavity, at the uterotubal junction, and 
inside the fallopian tubes (mucosal epithelium) where cells are stored and undergo chemical and mechanical changes 
17 1 . The cumulus oophorus, thin matrix covering the outside of the ovum, is also a viscoelastic (actin-based) gel 
181 ] . These complex fluids act as mechanical barriers for the selection of healthy spermatozoa and, together, lead to 
a decrease of six orders of magnitude in the size of the population of swimming cells, from tens of millions reaching 
the cervix to only tens approaching the ovum. 

Early rheological studies showed that the cervical mucus viscosity is about 100 Pas, with relaxation times on 
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111 ]. Subsequent measurements reported smaller relaxation times, between 1 and 10 s 
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24|]. These number vary during the menstrual cyclej^^ and depend strongly on the hydration of 
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24j and viscosities in the range 0.1 — 10 



the mucus [l{J. Given that the flagella of spermatozoa typically beat at a frequency between 20 and 50 Hz 0], 
spermatozoa swimming in mucus occurs therefore also at large Deborah numbers. The rheology of cumulus oophorus 
was characterized as well, with relaxation time of less than half a second for small deformation 18[ . The microstructure 
of cervical mucus was later examined by microscopy [26], and revealed a crosslinked gel with a typical length scale 
between polymeric filaments below one micron. 

The viscoelastic forces resulting from a complex mucus rheology significantly influence the waveforms, structure 
and swimming path of spermatozoa [271 ] . Early experimental study of the motion of human sperm cells in cervical 
mucus showed that, compared to motion in semen, their flagella beat with higher frequency, smaller amplitude and 
smaller wavelength, resulting in unchanged swimming velocity but straighter paths 28(. A subsequent study of the 
motion of bull spermatozoa in cervical mucus also reported lower swimming speed per unit beat than in a viscous 
fluid, as well as important wall effects 29], 

A last important issue regarding the motion of spermatozoa in complex fluids is called "hyperactivation" . It is a 
change in their flagellar beat pattern which occurs during their stay on the mucus-filled epithelium of the Fallopian 
tubes. The beat becomes asymmetric and results in higher force generation, helping the cells detach from the mucosal 
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epithelium, enhan cing their mobility in viscoelastic media, and facilitating penetration of the matrices surrounding the 
oocyte 30|, |31[ |32| . |33j . Experimental study of the motion of active and hyperactive spermatozoa in polymer solutions 



showed that hyperactive sperm swim faster than active (normal) cells as the elasticity of the solution increases (the 
viscosity remaining constant) and that hyperactive sperms swim in straight lines only in viscoelastic media [34 1. 



B. Modeling approaches in flagellar and ciliary fluid mechanics 



We now turn to the modeling of locomotion and transport in biological fluids. Much work has been devoted to the 
fluid mechanics of flagella and cilia in Newtonian flows, and the reader is referred to classical texts for reviews 0, 0, 
4, 35]. Quantitative approaches for ciliary flow include infinite models, where the envelope of the cilia tips is replaced 
by an infinite waving two-dimensional or three-dimension surface ll , 13611 , envelope models where an impenetrable 
oscillating surface is embedded on a finite three-dimensional object [3j],|38j], and sublayer models which detail the fluid 
forces at the individual flagella level 39, 40, 41 , 42 , 4^| . Classical modeling approaches in single-flagella hydrodynamics 
include resistive-force theory, where the fluid forces are described by a local anisotropic drag tensor and 



slender-body hydrodynamics [2j, 45, 46, 47, 48, 49, 50 1 where long and thin filaments are represented by a line 
distribution of flow singularities (see also Ref. 27j for an introduction to recent numerical work) . 

A number of modeling approaches for mucus transport and sperm locomotion in complex fluids have been proposed 
in the past. Respiratory mucus transport was quantified by two- fluid models, all showing that the mechanical proper- 
ties of the mucous layer are critical for efficient transport 
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55| . Taylor's original calculation of propulsion 
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by a waving sheet [jj, [56J was adapted to a variety of configurations relevant to a complex fluids, in particular in the 
case of interaction with oscillating walls 57|, interaction with a phenomenologically active material 1581 1. motion in 
a confined environment 
fluids 



6l| (see also related work on presitaltic transport of complex 
6 A l65lp . Numerical computations for the motion of flagellated cells interacting with deformable walls 
were reported [661 ] . showing swimming velocities decreasing with an increase in the wall elasticity. 

Only a few studies incorporated the viscoelastic properties of the fluid in a continuum-mechanics approach. The 
first effect of fluid elasticity was quantified using a retarded-motion expansion and assuming a second-order fluid 



671 ]. A more rigorous integral constitutive relationship was also used [681 ] . For zero Reynolds number, the swimming 
velocities reduce to the Newtonian velocities in both studies 67, fig}. A phenomenological modification to the equations 



of the resistive-force theory of slender-body hydrodynamics was proposed as a model for flagella motion in a linear 
viscoelastic fluid jf^. In that limit, the swimming velocity is unchanged from the Newtonian value, but the work 
done by the swimmer always decreases, and a complex fluid is seen to increase the swimming efficiency. A somewhat 
similar model was recently used to model changes in flagella waveforms in beating elastic filaments 7C | . 



C. Purpose of the study and structure of the paper 

In this paper, we propose to take a rigorous modeling approach to the fluid rheology in order to obtain quantitative 
information on transport and locomotion in a viscoelastic fluid in the biologically-relevant limit of large Deborah 
numbers. We consider Taylor's waving sheet as a idealized geometrical model of both spermatozoa locomotion (free 
sheet, quiescent flow in the far field) and ciliary transport (fixed sheet, finite induced velocity in the far field) in 
complex fluids. Although this is clearly a simplified model, it allows us to treat rigorously the effect of viscoelastic 
stresses on transport processes, and is therefore likely to shed some light on the biological examples detailed above. 

Using two models derived from the kinetic theory of polymers, we demonstrate that the viscoelastic nature of the 
fluid modifies in a non-trivial manner the kinematics and energetics of swimming of the sheet and transport of the 
surrounding fluid. Compared to the case of a purely Newtonian solvent, swimming and transport with the same gait 
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FIG. 1: Schematic representation and notation for the swimming sheet. 



always require more work and lead to lower speeds. Our results are therefore in stark contrast with those of Ref. 69 1 
where a linear viscoelastic material was considered. This arises because, generically, the net magnitude of locomotion 
and transport scale quadratically with the amplitude of the local oscillatory motion, and therefore nonlinear terms 
in non-Newtonian constitutive relationships can never be neglected. Surprisingly, we also show that our results are 
unchanged if wc consider other classical empirical constitutive relationships to describe the fluid dynamics. One 
implication of our results is that, although non-Newtonian fluids lead to a deterioration of the transport performance, 
they allow on the other hand biological systems to tune passively transport kinematics by modulating material 
properties. This feature is notably absent in Newtonian fluids and might suggest an advantageous biological control 
mechanism. 

The structure of the paper is the following. In Sjl] we present the setup and notation for the swimming sheet 
model. The swimming velocity of the free sheet (or induced far-field velocity for a fixed sheet) and its rate of work are 
then calculated in £ 11111 and §IVI for the Oldroyd-B and FENE-P models with a single relaxation time respectively. A 
discussion of the results and their relevance to biology is presented in Jvl In Appendix [Aj we re-derive the results for 
a second-order fluid. We then show in Appendix [51 and lO that other classical rheological models (Johnson-Segalman- 
Oldroyd and Giesekus) lead to the same results as those obtained in the main section of this paper. Furthermore, 
we show in Appendix [D] that our results are unchanged if we consider a sheet with traveling waves of both normal 
and tangential motion. Finally, the main results of this paper are generalized in Appendix [E] to the case where the 
viscoelastic fluid possesses more than one relaxation time. 



II. SETUP 

We consider in this paper the simplest model for an oscillatory swimmer, analogous to that proposed by Taylor in 
his original paper [l| . An infinite inextensible sheet of negligible thickness is swimming steadily in an incompressible 
viscoelastic fluid (see Fig. [T] for notations). In the frame moving at the (yet unknown) swimming speed, the vertical 
position of the sheet, described by the function y(x,t), oscillates in time according y(x,t) — asiri(kx — tot). This 
profile is a traveling wave of amplitude a, wavelength A = 2ir/k and period T — 2ir/u along the positive x direction 
(w.l.o.g both k and lo are assumed to be positive) and is responsible for the swimming motion. If the sheet is not 
free to move, a velocity opposite to the swimming velocity will be induced in the far field, and therefore the same 
calculation provides models for both locomotion and transport in a complex fluid. 

We nondimcnsionalizc lengths by 1/fc, times by ui, and velocities by the traveling wave speed, c = oj/k, so that 
the dimcnsionlcss vertical position of the sheet is written y(x,t) — esin(a; — t), where we have defined e = ak. In 
this paper, we will assume that the amplitude of the sheet is much smaller than the wavelength and will derive the 



swimming results in the asymptotic limit e <C 1 90(. 
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A. Kinematics 

Since the sheet is two dimensional, we will look to solve the problem using a streamfunction, ip(x, y, t), so that the 
two components of the velocity field, u, are given by u x — dij)/dy, u y — —dip/dx, and the incomprcssibility condition, 
V • u = 0, is always satisfied. The boundary conditions for the streamfunction arise from conditions at infinity and 
on the waving sheet. Let us denote by —Ue x the steady swimming speed of the sheet, where we use a minus sign 
because we expect swimming to occur in the direction opposite to that of the traveling wave as in the Newtonian 
case. In the frame moving at the steady swimming speed, the velocity field far from the sheet asymptotes therefore 
to Ue x , leading to the boundary conditions 

On the swimming sheet, the no-slip boundary condition leads to 

Vtb\, . , ... = ecos(a; — t)e x . (2) 

We look for a regular perturbation expansion in powers of e for both the streamfunction and the swimming velocity, 
in the form 

i/> = etp! + e 2 i> 2 + U = eU 1 +€ 2 U 2 + (3) 

and because of the symmetry e — -> — e, we expect swimming to happen at best at order e 2 . Expanding the boundary 
conditions on the sheet, Eq. @, we find that the leading order problem satisfies 

V ^ 1 l(a ! ,o) = cos(x~t)e x . (4b) 
while the boundary conditions for the stream function at order e 2 are 

V ^ 2 l( K ,oo) = U 2 e V ( 5a ) 

(5b) 



V H X , 0) = - sin ^-*) v (^ 



B. Kinetics 

We assume the Reynolds number for the fluid motion to be very small, and write therefore mechanical equilibrium 
of the fluid as 

Vp = V • t (6) 

where p is the pressure and t the (deviatoric) viscoelastic stress tensor. The pressure can be eliminated from the 
mechanical equilibrium by taking the curl of Eq. © and considering therefore the resulting vorticity equation, which 
is a fourth-order partial differential equation for ip, with boundary conditions given by Eqs. ^ and ([5]) for the first 
two orders in e. 

In order to close the system of equations and solve the pro blem we need to have a constitutive equation relating 
stresses to kinematics. Following classical monographs 

0,13, H uA l75l . l76l . 1771 178| | , we decompose the viscoelastic 



stress, r, into two components, r = t s + t p , where t s is the Newtonian contribution from the solvent, t s — i] s j, 
with 7 = Vu + Vu T , and t p is the polymeric stress contribution. Modeling choices govern the evolution equation for 
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III. MODERATE DEBORAH NUMBERS: OLDROYD-B 



As our first model, we consider non-Newtonian fluids governed by the Oldroyd-B constitutive equation with a single 
72, 73[ [zjl EE Izsl El, 79 1 . This is arguably the most famous constitutive equation for polymeric 



relaxation time [71 

fluids, and for two reasons. Firstly, it is very simple, and as so is very attractive. Secondly, it can be exactly derived 
from the kinetic theory of polymers in the case where the polymer molecules are modeled by elastic dumbbells, and 
therefore is physically sound. Oldroyd-B fluids agree reasonably well with experimental measurements of polymeric 
flows up to order one Weisenberg numbers (We = A7 ~ 1), but is known to possess a number of unphysical features 



at large values of extension rates when elastic dumbbells become infinitely elongated 7l|, [72 , 
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Model 



In the Oldroyd-B model, the polymeric stress, t p satisfies the upper-convected Maxwell equation 



r p + At* 



%1, 



(7) 



where A is the polymer relaxation time and r\ v the polymer contribution to the viscosity. The relaxation time measures 
the typical rate of decay of stress when the fluid is at rest and is zero for a Newtonian fluid. The upper-convected 
derivative, defined for the tensor A as 



dA 
~dt 



u • VA - (Vu T • A + A • Vu), 



(8) 



is the rate of change of A calculated while translating and deforming with the fluid. It is then straightforward to see 
that Eq. (O can be manipulated to obtain the Oldroyd-B constitutive equation for the total stress, x, as 



t + Ai t= r)[y + A 2 7] 



(9) 



where rj — rj s + rj p , \\ = A, and A2 = 77^ A/ 77 < Ai, and which is the form we will use to solve the swimming problem. 
Note that A2 is termed the retardation time scale and is equal to the rate of decay of residual rate of strain when the 
fluid is stress-free. 

We now nondimcnsionalizc stresses by tjuj and shear rates by ui, leading to the constitutive relationship 



t + Dei t= 7 + De 2 7, 



(10) 



where we have defined the two Deborah numbers Dei = Aiw and De 2 = A 2 w and where we use the same symbols for 
convenience. 

B. Order e solution and rate of work 



We look for a solution for the stress as a regular perturbation expansion in powers of e, 



T — €Ti + 6 T 2 



and the leading order solution of Eq. (|10[) is written as 



^ cVi . ^ 971 
n + Dei— = 7i + De 2 ^- 
ot at 



(11) 



(12) 
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Since Eq. (|12p is a linear equation, we can first calculate its divergence and then its curl to obtain the equation for 
the leading order vorticity, and we get 



1 + De 2 ^ ) VV, =0. 



(13) 



Given the boundary conditions from Eq. (U) , and using Fourier transforms, it is straightforward to see that ipi is the 
same as for the Newtonian case, that is 



ifo (x, y, t) = (l + y)e~ y sm(x -t), Ui=0. 



(14) 



At leading order in e, the flow field and swimming kinematics are therefore unchanged by the presence of viscoelastic 
stresses. The rate of work of the sheet, however, is modified. Let us use Fourier notations, and write 



-Vp-tx 



Using similar notations for the stress and rate-of-strain tensors, we obtain the first-order stress as given by 



1 + «De, 



Tl 



-7r 



(15) 



(16) 



1 + iDei 

The leading-order rate of viscous dissipation, written as W = e 2 W 2 + is equal to the volume integral of t\ : ■y 1 on 
both sides of the sheet, which average value in Fourier notation is the volume integral of 



[Tl ■ 7i) = ^i T i = 7il = 77 ( — ^2 ) [7i : 7iL 



(17) 



2 * llJ 2 V 1 + Dc 2 

where (...) denotes the average over one wavelength of the sheet. As a consequence, the leading order value of the 
rate of viscous dissipation, equal to the rate of work done by the sheet, is given by, in dimensional variables and per 
unit length in the z direction (see Fig. [J), 



W 2 _ 7? 1 + DciDc 2 
Wn ~~ t] N l + De? ' 

where WV denotes the leading-order rate of working in a Newtonian fluid of viscosity tjn. 



(18) 



C. Order e solution and swimming velocity 



We now turn to the determination of the swimming velocity of the sheet. The order e 2 component of Eq. (|10p is 
given by 



1 + Dex^J T 2 -(l + Dea^J 7 2 = D Cl (Vuf • n + n • Vu x - Ul • Vti) 

-De 2 (Vuf • 7i + 7i • Vux - u x • V7O 



(19) 



The only relevant part of the streamfunction ^2 is its mean value, as this is the only one possibly contributing to 
a non-zero swimming speed in the far field. Using Fourier notations and Eq. I|16p . we get the mean value of the 
constitutive equation, Eq. (|19p . as given by 

(20) 



The right-hand side of Eq. (|2"U)) can be calculated using the leading-order kinematics, given by Eq. (fT5"|) . We find 



Ui = 



-ly 
■0-+V) 



-Ve—tX 



, VUl = 



-y + 
i(y-i) y 



7i = 



-2y 2iy 
2iy 2y 



(21) 
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leading to 



6y 2 -2y-l Dei(l + 2j/ - 2y 2 ) 
D ei (l + 2y-2y 2 ) 2j/ 2 + 2y + 1 



(22) 



We then take the divergence of Eq. (f2"2"|) , and the curl to eliminate the pressure, and obtain the equation for the 
strcamfunction 

4 D ei (De 2 -Dei) d 2 

1 + Dc 2 dy 2 

with boundary conditions given by the time-average of Eq. ([5]) 

d{ih) 



[e- 2 v(l + 2y - 2y 2 )} , 



dy 
dy 



(oo) - C/ 2 , 

(0) - i 



The solution to Eq. (|23[> is, using the boundary condition given by Eq. (l24b|) . 



so that the swimming velocity is given by 



U N 



1 + D ei De 2 
1 + Del 



(23) 

(24a) 
(24b) 

(25) 

(26) 



where Un = 1/2 is the leading-order swimming velocity obtained in a Newtonian fluid. The velocity of the free 
sheet (or, alternatively, the flow speed induced in the far-field for a fixed sheet) is therefore always smaller than the 
Newtonian value (U < Un). 



IV. HIGH DEBORAH NUMBERS: FENE-P 



As was discussed above, the Oldroyd-B constitutive equation is a very attractive model for polymeric fluids, but it 
is not valid for large extension rates. In fact, in a purely extensional flow, its extensional viscosity blows up at a finite 
value of the extension rate. For the setup considered in this paper, the flow near the sheet is a combination of shear 
and extensional flow, and we are therefore not allowed to use the results obtained in the previous section, Eqs. fl8]) 
and (|26|) . in the biologically-relevant limit of large Deborah numbers. In order to remedy this issue, we consider in this 
section a more advanced constitutive relationship, the FENE-P model 
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74 



75 



77 



of the FENE-P model is very similar to Oldroyd-B, with the difference that the entropic sp ring s modelin, 
molecules have finite extension and can therefore sustain arbitrarily large strain rates 
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80j . The physical picture 
olymer 
1. The 
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model can be exactly derived for a solution of finite-extension dumbbells expect for one empirical averaging step (the 



Peterlin approximation [81|, |82j, at the origin of the "P" in FENE-P). 



A. Model 



In the FENE-P model, the polymeric stress is given by 

t* = G[/(k)k-1] 



(27) 
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where G has units of stress, 1 is the identity tensor, and k is the dimensionless dumbbell covariance tensor, which 
satisfies the evolution equation 

/(k)k + Ak = 1, (28) 
where A is the (single) polymer relaxation time. The function / is defined as 

/(«) = — W ( 2g ) 

where L is a (large) dimensionless number measuring the maximum relative extension of the polymer molecules. In 
the limit where L — ► oo, this model reduces to the Oldroyd-B model, Eq. ([7|), with r\ v = G\. 

As in the previous section, shear rates and stresses are nondimensionalized by oj and r\uj respectively, and the value 
of the total viscosity, 77, will be defined below. All symbols will now refer to dimensionless variables. 

B. Expansion 

At equilibrium, we obtain the order e° solution, with no flow (uo = 0) and 

160 = 7-^-3-1. /(*>) = /o = 1 + t p = Q. (30) 
We then look for a regular perturbation expansion in the form 

K - Kq = 6Kl + e 2 K 2 + (31) 

and therefore 

/(«)-/o = e/i+e 2 /2- (32) 

Using Eq. (|3"T|) , we see that 

, (3 + L 2 )\ , (3 + L 2 )\ , (3 + L 2 )\ 

h = jje tr («i)! /a = 7^ tr(Ka) + — 775 — tr(«i). (33) 

C. Order e solution and rate of work 

At order e, we have 

/i«o + /o«i + Xu} ~^f = Aw [Vuf • /«o + «o • Vui - ui • Vk ], (34) 

and the last term on the right-hand side of Eq. (fM)) is zero as k q is constant. Using the fact that k = l//o, Eq. 
becomes 

fiKo + /o«i +Aw^i = ^7 X , (35) 
cr Jo 

and by taking the trace of Eq. ((35|) we find tr«i = /1 = 0. It is then natural to define the first Deborah number as 

Dex = ^, (36) 
Jo 

and since (/o - 1) « 1, this definition is essentially equivalent to the definition Dei — \lo for the Oldroyd-B model. 
We then obtain, using Fourier space notations, 

1 Dei t 
fo 1 + zDei 
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The equation for the dimensionless polymeric stress is jjwtJ = G/0K1, so that the total stress, n = t\ + t\ , is given 
in Fourier space by 

TjUTx = r ]s ujj 1 + G i r (38) 

1 + iDci 

In the limit of small Deborah number Dei, the expression above allows us to define without ambiguity the polymer 
viscosity, and therefore the total viscosity and second Deborah number, which are given by 

GX r/ s 
Vp = -*-) V = Vs + rjp, De 2 = Dei — - (39) 
fo V 



The first-order non-dimensional stress is finally by 



1 + iDe 2 ~ , , 

^ = T~MDef 71 ' ^ 



This equation is identical to the one obtained for the Oldroyd-B model, Eq. (fl6|) , and therefore we find that, provided 
that the Deborah number and polymer viscosity are defined using Eqs. (|36[) and (|39[) . the leading-order flow kinematics 
as well as the rate of work of the sheet are the same as for the Newtonian case, and given by Eqs. (p~5|) and ([18]) . 

D. Order e 2 solution and swimming velocity 

At next order, the evolution equation for n becomes 

/2«o + /o«2 + fo^i-^f = /oDei[Vu^ ■ kq + k - Vu 2 + Vuf • ki + «i • Vui - m ■ V«i], (41) 



and the polymer stress is given by 



T 



p _ ^ Dei - De 2 

2 — t-.„2 



V De? 



[/ 2 «o + /o«a] , (42) 



where have used that G/r/uj = (Dei — De 2 )/De^. Combining Eqs. (|4"Tj) and (14"2")) . and exploiting the fact that Kq = l//o, 
we obtain the expression for the average of the second-order stress as 

De 2 \ . / De 2 \ T 



<r^) = ( 1 - (7 2 ) + /o ( 1 - (Vu? • Ki + Ki ■ Vui - Ul • Vki), (43) 



and therefore, with r 2 = i" 2 + f 2 , we obtain 



(t 2 > - <7 2 ) = /o ^1 - (Vu[ • ki + ki • Vui - ui • Vki). (44) 

Using Fourier notations we find 

<T 2 > - (7 2 ) = H { y (l - §|) (Vui* ■«!+«!■ V5J - uj • Vki) j , (45) 
so that, exploiting the result of Eq. ([37)1 . we obtain 



(t 2 > - <7 2 > = » { 2^^) ( V ^* ' 71 + ^ ' - 3 * ' } ' (46) 

This expression is identical to the one obtained for the Oldroyd-B model, Eq. (|20|) . and therefore the swimming speed 
of the sheet is the same as well, given by Eq. 
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V. DISCUSSION AND IMPLICATIONS FOR BIOLOGICAL TRANSPORT 

In this paper, we have shown that the presence of a non-Newtonian fluid modifies significantly the kinematics and 
energetics of transport and locomotion by a waving sheet as compared to that in a Newtonian fluid. Specifically, we 
have demonstrated that the leading-order swimming velocity of a free sheet (or, alternatively, the induced far-field 
velocity field when the sheet is fixed) is given by 

U _ 1 + Be 2 r] s /r] 
Thi ' 1 + De 2 ' 

where 

1 



(47) 



U N = ^a 2 k, (48) 

is the leading-order (dimensional) velocity obtained in a Newtonian fluid; next order terms in Eq. (|48p are 0(u>a^k 3 ). 
With the appropriate definition of the Deborah number and the polymer viscosity, this result was shown to be valid 
for two constitutive relationships, Oldroyd-B ( £11111) and FENE-P f ^IVj) . but other models such as Johnson-Segalman- 
Oldroyd (see Appendix IB"]) and Giesekus (see Appendix [C|l lead the same result for the swimming speed. Furthermore, 
as shown in Appendix [D] the result is also unchanged if the wave of displacement propagating along the sheet includes 
not only normal but also tangential motion. The generalization of Eq. (|47|) to the case of many relaxation times is 
presented in Appendix [El 

Since the zero-shear-rate viscosity of the polymeric fluid, T) = rj a +?]p> 1S larger than the viscosity of the Newtonian 
solvent, T] s (and in most cases, much larger, see JQ): we always have U < Um- Locomotion with a given gait is 
therefore always slower in a non-Newtonian fluid. In the biologically-relevant limit of large Deborah numbers, the 
velocity asymptotes to a constant value, equal to Un^s/Vi and therefore the swimming speed is significantly reduced 
from the Newtonian value. The dimensional dependance on the parameters of the swimming gait, Eq. (|48|). suggests 
that such reduction in velocity could be compensated by an increase in the beat frequency (still consistent with the 
De ~3> 1 limit) and a decrease in the wavelength, both of which are observed for spermatozoa motion in cervical mucus 
. An increase in the beat amplitude would also increase the swimming velocity, however sperm cells are observed 



to swim in cervical mucus with decreased amplitude [281 ] . 

The leading-order rate of work done by the waving sheet, as compared to the work Wjv in a Newtonian fluid of 
viscosity 77 at, is given by 

W = n l + De^/r? 
Wn W 1 + De 2 

Notably, this result is identical to that given by forced oscillations of the Jeffreys spring-dashpot model: a dashpot 
of viscosity 7] s in a parallel with a Maxwell model, i.e. of a dashpot of viscosity rj p in series with a spring of constant 
k p = i] p /X. 

It is reasonable to consider two different cases when discussing Eq. (|49|) . Firstly, we can compare the motion in the 
polymeric fluid with the motion in the solvent alone, which requires a rate of work W s . In that case, we need to pick 
Vn = Vs, and Eq. ([49]) becomes 

W y/ Vs + De 2 

W. = "T+De^' (50) 

Swimming in a polymeric fluid requires therefore always more work than swimming in the solvent alone (W > W s ); 
note however that in the limit of large Deborah numbers, we get W « W s . We can also compare the rate of viscous 
dissipation with that occurring in the Newtonian fluid which has the same zero-shear-rate viscosity as the polymeric 
fluid, with rate of work W r) . In that case, we choose t/n = rj and obtain 

W _ 1 + DeV/?? 

W v ~ 1+De 2 ' [b) 
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As a difference, we see here that compared to a Newtonian fluid with the same viscosity, it is energetically advantageous 
to swim in a fluid that has some elasticity (W < W v ). 

The mechanical efficiency of the swimming motion, defined as the ratio of useful work to total work, is given by 

£~^r> (52) 

where the units in Eq. (|52p are correct as W is a rate of work per unit length in the spanwise direction (z direction, 
see Fig. [1]). Given Eqs. (gT]) and (gSJ, we obtain 

e = i + DeV, 

£n l + De 2 

when £n is the efficiency for the swimming motion in the Newtonian fluid. Again, we obtain £ < £n and swimming 
in a non-Newtonian fluid deteriorates the mechanical efficiency. 

At first, these observations seem to suggest that there are only disadvantages to swimming in a complex fluid, as 
compared to a Newtonian fluid. However, we would like to argue here that non-Newtonian fluids could be exploited 
by biological systems to their advantage because viscoelastic stresses allow a tuning of the kinematics of transport 
and locomotion in a manner which is not possible with a Newtonian fluid. Indeed, in the case of a Newtonian fluid, 
the swimming velocity of the sheet is given by Eq. (|48|) . and it is only a function of the swimming gait parameters. 
The only way to modify, for example, the speed at which cilia transport a Newtonian mucus would be to change 
the amplitude, frequency or wavelength of the ciliary waves, and a viscosity change would only impact the necessary 
rate of work. However, as Eq. (|47ll shows, the transport kinematics in a viscoelastic fluid are also controlled by the 
mechanical properties of the fluid, namely its viscosity and relaxation time. This new feature allows the transport to 
be tuned without having to modify the parameters of the swimming gait, a feature which is notably not present when 
the fluid is Newtonian and could potentially be important. 

This observation suggests the possibility, for a biological system, to modulate transport and motility passively. For 
example, tuning the rheology of respiratory mucus would influence the rate at which it is transported and evacuated 
without having to control the cilia themselves. Similarly, varying the rheological properties of the cervical mucus would 
allow the female reproductive tract to select appropriately motile spermatozoa without requiring a modification of 
the waveforms and beat frequency of their flagella. A number of experimental observations support this idea. The 
rheological properties of cervical mucus vary during the menstrual cycle, and measured "sperm penetrability" peaks 
during the ovulatory phase when the viscosity of the mucus is at a minimum 



25. 



8: J >. 



This is consistent with our 

result dU/drip < 0, easily inferred from Eq. (|47p. The elasticity of the fluid has also been reported to be an important 
factor 



85, 



86, 



87 1 and transport usually decreases with an increase in the elasticity [lj, |23[ , an observation which is 



consistent with our result dU/dDe < (see Eq. |4"T1) 91]. Similar results were obtained in numerical simulation of 

ri . ri 

mucociliary transport [511 ] and spermatozoa motion in deformable domains [66] . 

In conclusion, we have revisited in this paper Taylor's swimming sheet calculation in the case where the fluid is 
non-Newtonian. Although it is clearly a very simplified model, the results obtained in this paper could be extended to 
more complicated geometries and wave profiles, and could therefore provide a first step towards a rigorous quantitative 
approach to the energetics and kinematics of locomotion and transport in complex biological fluids. 



APPENDIX A: RETARDED-MOTION EXPANSION 



In this section, we re-derive the result first obtained in Ref. [671 ] and show that a waving sheet in a second-order fluid 
swims at the same speed as in a Newtonian fluid. We also show that the rate of working of the sheet is unchanged. 

A second-order fluid is the first term in a systematic asymptotic expansion of the relationship between the stress 
and the rate-of-strain tensor for small Deborah numbers and is valid in the limit of slow and slowly varying flows 



0,B 



73 



75 



76 



77] . The instantaneous relationship is given by 



1 v 

t = voi - 2* 1 ^ +* 2 v7 ■ 'V)' 
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(Al) 



where and ^ are the first and second normal stress coefficients. We nondimensionalize stresses by r/ouj, shear 
rates by u>, and the relationship becomes 



r = 7 - Dei 7 -De 2 (7 • 7), 



(A2) 



where Dei = ^iu>/2rjo and De 2 = — Vt^w/^o are the two characteristic Deborah numbers (both are expected to be 
positive). 

At first order, Eq. (|A2|) becomes 



Tl 



and therefore the strcamfunction satisfies 



1-Dei|bv 



d 



1-Dci- )V*Vi=0. 



(A3) 



(A4) 



Similarly to the result of i jlll B[ we see that the solution is the same as that obtained by Taylor for motion in a 
Newtonian fluid as given by Eq. (|T4"]) . The first-order component of the stress tensor is given by 



Ti = (1 - iDei)7i, 
and therefore the rate of work is equal to the volume integral of 

{ti ■ 7i) = 2 K { ?1 : 7*} = 7i ■ 7i, 

and is the same as for the case of a Newtonian fluid (with equal viscosity 770). 
At second order, Eq. (|A2j) is given by 

r 2 - (l - De iJ^ 72 = Dci(Vuf -7i +7i ' Vui -ui • V71) - De 2 (7i -7^, 

and its mean value is given by, using Fourier notations, 

(t 2 ) - (7 2 > = 1® {Dei (Vuf* ■ \ + \ ■ Wul - u{ • V^) - Ve 2 fii ■ . 
Using Eq. (fTS"]) , we get 



(A5) 



(A6) 



<T2) - (72) = e 



-2y 



(6y 2 - 2y - l)Dei - 4y 2 De 2 

(2y 2 + 2y + l)De! - 4y 2 De 2 



for which we can take the divergence and then the curl to obtain 

V 4 (^ 2 ) =0. 



(A7) 



(A8) 



(A9) 



(A10) 



The equation and the boundary conditions are the same as for the Newtonian problem, and therefore the swimming 
velocity is unchanged. 
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APPENDIX B: JOHNSON-SEGALMAN-OLDROYD 

The Johnson-Scgalman-Oldroyd constitutive relationship is an empirical model capturing non-affine motions in 



the polymeric fluid 88J . It is formally similar to the Oldroyd-B constitutive relationship but uses the (more gen- 



eral) Gordon-Schwalter convected derivative. It includes other simple models such as Oldroyd-B, Oldroyd-A, and 
corotational Maxwell as special cases 

HHH- 



The constitutive relationship is given by 

□ D 

T + AiT=f7[7 + A 2 7], ( B1 ) 
where the Gordon-Schwalter convected derivative is 

-r=T+-(7--r + T-7), (B2) 

and where a is a dimensionless parameter related to second normal stress differences and leading to shear-thinning 
behavior in steady shear. Using the same non-dimensionalization as for Oldroyd-B, we obtain the dimensionless 
constitutive relationship 

□ D 

r + Dei r= 7 + De 2 7 . (B3) 

Since the linear part of Eq. (|B3|) is the same as for Oldroyd-B, it is straightforward to see that the problem at order 
e is also the same as for Oldroyd-B, so the rate of working of the sheet is given by Eq. (|49p. The change in the 
convected derivative will affect at most the second-order solution. The constitutive equation at second order has two 
terms added as a result of the change in the convected derivatives 

I + Dei^ t 2 - (l + Vc 2 ^) 7 2 = De x (Virf • t x + n • Vuj - u t • Vt x ) (B4) 

-^ De i(7i -Ti + ri -7!) 
-De 2 (Vuf • j 1 + 7 X • Vui - ui • V7J 
+aDe 2 (7! -71), 



dt V at 



which leads to the modified equation for the mean stress 

(T2>-(72) = 



a ( De 2 — Dei 



4V i+Dzr'^^-™ 

The new term on the right-hand side of Eq. (|B5[) is given by 

7i -7i = 16y 2 e~ 2y l. (B6) 

The tensor given in Eq. (|B6p is proportional to the identity tensor, 1, and can therefore be absorbed in the definition 
of the fluid pressure. Consequently, this new term does not modify the equation for the streamfunction as given by 
Eq. ([23| and the swimming velocity is the same as for Oldroyd-B (Eq. |47|) . All Oldroyd-like models lead therefore to 
the same results. 



APPENDIX C: GIESEKUS 



As our last model, we consider the Giesekus constitutive relationship 
using molecular ideas for polymer networks 

0, H H H H 



7G 



, which, although empirical, was derived 



77] and is therefore directly relevant to complex 
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biological fluids. In this model, the polymeric contribution to the stress satisfies the evolution equation 



v A 
t p + Xt p + a—(r 



p . 



(CI) 



where a is a small dimensionless positive number related to the second normal stress coefficient in steady shear. The 
equation for the total stress r is then found to be 



r + Aj 



t — a — (7 • t + t ■ 7) 

Vp 



a— (t • r) 



7 + A 2 



a— (7 -7) 

Vp 



(C2) 



where, again, -q = rj s + rj p and Ai = A and A 2 = Xrjs/r]. Compared to the Oldroyd-B constitutive equation, three 
new quadratic terms appear. The new terms of the from (7 ■ 7) and (7 • t) are formally similar to the extra terms 
appearing in the Johnson-Segalman-Oldroyd constitutive relationship, Eqs. ()B1|) and (|B2|) . and do not modify the 
swimming velocity (see Appendix [B| . We therefore only need to quantify the influence of the term quadratic in the 
stress, last term on the left-hand side of Eq. (IC2jl . In non dimensional and complex notations, we still have at first 
order 



1 



iDe 2 ~ 
7i> 



1 + zDci 



(C3) 



so that the expression for the rate of work is unchanged, and the steady contribution to the second-order streamfunction 
arises from the new term 



(t% -Ti) = |(tx -tI+tI -n) 



1 (1 

4 



De, 



1 + Det 



(7i -7i +7i -7i] 



(C4) 



As for the Johnson-Segalman-Oldroyd model, this term does not modify the equation for the mean value of the order- 
two streamfunction, and therefore does not impact the swimming kinematics. As a consequence, our final results, 
Eqs. (|4T|) and (|4*5|) , are also valid for the Giesekus model. 



APPENDIX D: ALLOWING TANGENTIAL MOTION OF THE SHEET 

Here, we show that the main results of this paper, Eqs. (|4"T)) and (|4"5)) , are unchanged is we allow the wave of 
displacement along the sheet to include not only normal but also tangential motion. In the reference frame moving 
with the sheet, we write the dimensionless position of material points, (x m ,y m ) on the sheet as 

x m = x + eacos(x — t — (ft), y m — ebsm(x — i), (Dl) 

where a and b are dimensionless and 4> is the phase difference between the tangential and normal motion. The 
boundary condition in the far-field is unchanged, 

V ^,oo) =f7e v> (° 2 ) 
whereas the boundary conditions on the sheet are modified and given by 

V^L n = e(bcos(x — t)e x + asin(x — t — <ft)e y ). (D3) 

At order e, they become 

Uiey, (D4a) 

b cos(x — t)e x + a sin(x — t — (f>)e y , (D4b) 



V H*,°°) = 



16 



and at order e 2 we obtain 



(x,0) 



U 2 e v , 

—a cos(x — t — 0)V 



dx 



bsin{x -t)V 



(x,o) V dy 

For the Oldroyd-B model, the problem at order e is straightforward to solve and we obtain 

& = i[6 + (6 + ae^e^e"^ = ib[l + (1 + c)?/^^", 



(x,0) 



(D5a) 
(D5b) 

(D6) 



where we have defined, c = ae^/fe. Since the equation quantifying the mean value of the rate of work, Eq. (|17p . is 
not modified, the final formula for the rate of work of the sheet is unchanged, and given by Eq. ([49]) . We can then 
compute 



Ui = b 



ic — i(l + cjy 
-\-{l + c)y 



(D7) 



Vui = b 



c — (1 + c)y i + + c)y 

+ 2c) + i(l + c)y -c+(l + c)y 



(D8) 



Vu* = b 



c*-(l + c*)y -i-i(l + c*)y 
i(l + 2c*) - i{\ + c*)y -c* + (1 + c*)y 



(D9) 



and 



T* 



c*-(l + c*)y i(l + 2c*)-i(l + c*)y 
-i-i(l + c*)y -c* + (l + c*)y 



e- y e lx , 



ii=b 

Using the previously-derived relationship 



(r 2 > - <7 2 > = 
we obtain, after some algebra, 



2c- 2(1 + c)y -2ic+ 2i(l + c)y 
2ic+ 2i(l + c)y -2c + 2(l + c)y 



Dei - De 



2(l + zDe 



7i + 7i • VuJ - u{ • V-Vi 



V 4 ( ^ 2) = Dei (De 2 ^ Dei ) ^ 2 ^_ ^- 2j/(1 _ 3cc * + 2y(1 + 2c + 2c * + 3cc * } _ 2y2(1 + c)(1 + c *))] 



associated with the boundary conditions 



dy 



(oo) = t/ 2 , 
1 



(0) = -(b 2 + 2abcos<t>~ a 2 ), 
ay 2 



Eq. (|D13|) is easily integrated three times and we obtain 



U 2 = \{b 2 + 2a6cos</> - a 2 ) + l De i( De 2 De i) 6 2 (1 + c + c *_ cc * }j 
2 2 1 + De! 



(D10) 
(Dll) 
(D12) 
(D13) 

(D14a) 
(D14b) 

(D15) 
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(D16) 



and since b 2 (l + c + c* — cc*) = (b 2 + 2a6cos0 — a 2 ), we get 

U 2 _ 1 + DeiDe 2 
U N ~ 1 + Del ' 

where t/jv is the swimming velocity in a Newtonian fluid, 

U N = ^(b 2 + 2abcos<j>-a 2 ). (D17) 

Allowing tangential motion does therefore reduce the swimming speed by the same amount, and the final results of 
the papers, Eq. (f47|) and (|49|) . remain also valid in this case. 



APPENDIX E: MULTIMODE OLDROYD B MODEL 

In this section, we generalize the results of the paper to the case where the constitutive relationship for the 
viscoelastic fluid possesses more than one relaxation times [711.172, In that case, the polymeric contribution 



to the stress tensor is given by 

N 

r P = Y, T r ( E1 ) 
i=i 

where each mode satisfies an upper-convected Maxwell equation 

rf + A^ = % -7, (E2) 

characterizes by a relaxation time Xj and viscosity rjj. At low shear rates, the fluid has a steady-shear viscosity given 
by 

JV 

V = Vs+Y2 1 h- ( E3 ) 
i=i 

We now nondimcnsionalize stresses by rju, and shear rates by uj, and solve for the flow order by order. Using Fourier 
transforms, we obtain, in dimcnsionless variables, 



1 



where De^ = XjW. The total stress becomes therefore 



and the leading-order kinematics are unchanged from that obtained in the case of a Newtonian fluid, as given by 
Eq. (fT5|) . The rate of work done by the sheet is equal to the volume integral of 

' 1 r\ •f^ 77 1 + zDej 1 



and therefore we obtain 



W N Vn 
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where, as above, WW is the rate of work in a Newtonian fluid of viscosity t/n- In the case of a single relaxation time 
scale, Eq. (|E7|) is the same as Eq. ([4"9")) with rjx = T] p = r) — T] s and Dei = De. 

At next order in e, the mean value of the dimensionless form of Eq JE2| leads to 

<t? 2 > - ^<7 2 ) = De^Vu? • rl, + rl x ■ Vm - Ul • Vr^), (E8) 



which becomes, in Fourier space, and using Eq. (|E4 



(E9) 



We then sum up Eq. (|E9|I for each mode, and add the Newtonian contribution from the solvent, to obtain 

w -<*>-« ( (± 1^52-1 (var ■ if. + 5, • vs; - s; ■ v+, 



Using the matrices given in t^III C[ we obtain 

N 



(T2) - (7 2 > = e" 



I r] 1 + iDej 



6y 2 — 2y — 1 i(l + 2y - 2y 2 ) 
z(l + 2y-2 2 / 2 ) 2y 2 + 2y+l 



(E10) 



(Ell) 



and therefore 



(T 2 )-(7 2 )=e^E f T 



»» 1 + De" 



6y 2 -2y-l De^ (1 + 2y - 2y 2 ) 
Be J (l + 2y-2y 2 ) 2y 2 + 2y + l 



(E12) 



After taking the divergence and curl of Eq. (|E12[) , we obtain the equation for the mean streamfunction as 



V 4 (^ 2 ) = - 



N T7 De 2 



Def / dy 2 



[e^(l + 2y-2y 2 



(E13) 



and calculations similar to those presented in i jlll CI lead to the swimming speed of the sheet 

TT N T\ 2 



u N 



r-f V 1 



U 1 1 + De J 



(E14) 



In the case of a single relaxation time, the result of Eq. (|4T|) is recovered. Finally, the mechanical efficiency is given 
by 



£ 



N 2 



— = i - v — - 
£n f-f v i 



1 V l + Dej 



(E15) 
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